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Abstract 

We have carried out extensive equilibrium molecular dynamics (MD) simulations to study the 
structure and the interfacial properties in the liquid-vapor (LV) phase coexistence of partially 
miscible binary Lennard-Jones (LJ) mixtures. By analyzing the structural properties as a function 
of the miscibility parameter, a, we found that at relatively low temperatures the system separates 
forming a liquid A-liquid B interface in coexistence with the vapor phase. At higher temperatures 
and, < a < 0.5, we found a temperature range, T^(a) < T* < T^^^^^^a), where the liquid 
phases are wet by the vapor phase. Here, T^{a) represents the wetting transition temperature 
(WTT) and T^^^g{a) is the consolute temperature of the mixture. However, for 0.5 < a < 1, no 
wetting phenomenon occurs. For the particular value, a = 0.25, we analyzed quantitatively the 
T* versus p* , and P* versus T* phase diagrams and found, T* ~ 1.25, and T^ons — 1-25. We 
also studied quantitatively, as a function of temperature, the surface tension and the adsorption 
of molecules at the liquid-liquid interface. It was found that the adsorption shows a jump from a 
finite negative value up to minus infinity, when the vapor wets the liquid phases, suggesting that 
the wetting transition (WT) is of first order. The calculated phase diagram together with the 
wetting phenomenon strongly suggest the existence of a tricritical point. These results agree well 
with some experiments carried out in fluid binary mixtures. 



PACS numbers: 68.03.Cd 68.03.Hj 68.05. Cf 68.08.Bc 



I. INTRODUCTION 



Wetting phenomena occurs very often in daily life and has a strong interdisciplinary 
character. It is of great relevance for fundamental areas of knowledge such as physics 2 1 , 
chemistry j^, biology ^, and several other applied sciences P| as well as technology 6|. 
The occurrence of wetting phenomenon is usually associated with the existence of three 
phases, at least one of which is liquid and no more than one phase is solid. In most practical 
situations the solid phase is wet by the liquid phase and the disfavored phase is the vapor. 
This is expressed in terms of surface free energies as, < 75^. Although the vapor phase 
is disfavored, a solid-vapor coexistence can be possible and in such case it is said that the 
liquid "partially wets" the solid. However, as the liquid phase is further adsorbed by the 
solid, it is possible that the contact of the vapor with the solid be excluded by the creation of 
a liquid layer between them. In such a situation one says that the liquid "completely wets" 
the solid. The surface free energies involved in the wetting process are related by Young's 
rule, jsv = ISL + 7Lycos^^, with 6 the contact angle. For partial wetting, < < 90°, and 
for complete wetting, 6' = 0. Usually the presence of a wall in most of the theoretical as well 
as experimental studies complicates a detailed and precise investigation of the interfacial 
structures. This is due to the fact that the atomic interactions between the wall and the 
fluid components play an important role in the wetting phenomena. Fortunately, binary 
liquid mixtures offer a great opportunity to investigate in detail wetting phenomena since 
they involve only fluid phases in coexistence. An important number of experimental studies 



mm 



of interfacial wetting have studied only a small subset of possible binary mixtures 
0,0|. Wetting and prewetting phenomena in which one of the phases is solid has also been 
studied by means of Monte Carlo simulations using a 3D lattice-gas model ^12]. A 12-6 LJ 
potential modeling argon and a 9-3 L J potential modeling a CO2 covered solid wall have been 
chosen to study the wetting transitions by means of Monte Carlo simulations ^|. However, 
to the best of our knowledge, there are not numerical simulations that explicitly study the 
interfacial behavior and the associated wetting phenomena in fluid mixtures in absence of a 
wall. In this paper we apply a well established methodology using MD simulations 
llfil . IitI Q to study interfacial properties and surface phenomena in a partially miscible fluid 
mixture. We consider a model binary LJ fluid mixture in which the attractive part of the A-B 
interactions is weaker than the A-A and B-B interactions. By studying the properties of the 
density profiles of the fluid phases as a function of temperature and miscibility parameter, 



we show clear evidence that the vapor phase spontaneously excludes the liquid A-liquid B 
interface at and above the wetting transition temperature, T^. That is, the vapor phase fully 
wets the liquid A-liquid B phases. This phenomenon occurs when the miscibility parameter 
a is in the range < a < 0.5 and in the temperature region T^{a) <T* < T^Q^g(a). Here, 
T* = {kBT)/e, is the reduced temperature, T*{a), is the wetting transition temperature 
(WTT), and Tconsl*^)' consolute temperature of the mixture. By analyzing the 

structural properties of the mixture as a function of a we estimate the wetting phase diagram 
T^(a). A further quantitative analysis, for a = 0.25, of the surface free energy 7(T*) and 
the adsorption, r(T*), at the LL interface indicates that the wetting transition is of first 
order. 

The layout of the remainder of this paper is as follows. In section |n] we introduce 
the model of the symmetric binary mixture, in section IIIII we explain the details of the 
simulations, then in section IIVI we present and discuss the representative results of the 
extensive MD simulations. Finally, we end with the conclusions in section 

II. THE MODEL 

The model binary mixture studied in this paper consists of fluids A and B made up of 
spherical molecules of the same size, a^A = <^bbi and at concentrations of 50% each. The 
interaction between molecules of the same type is through the classical 12-6 LJ potential. 
However, the interaction between molecules of fluids A and B is given by the LJ potential. 



with tAA = ^BBi and uab the parameter that controls the miscibility of the two fluids. For 
the sake of simplicity from now on we will use the short notation, uab = ol. Notice that 
when, a = 0, we obtain two independent single LJ fluids while in the opposite case, a = 1, 
the system reduces to a single LJ fluid. By choosing, < a < 1, the attractive part of 
the A-B interactions becomes weaker than that of the AA and BB interactions, and then. 




(1) 



where the mixing rule is defined by 



0"AB = -^{(^AA + (^Bb), 



eAB = OiAB^AA, 



(2) 



the liquid phases are immiscible in a wide range of temperatures. Thus, one can obtain the 
coexistence of three fluid phases: liquid A-liquid B and the vapor. 

III. DETAILS OF THE SIMULATIONS 

We have carried out extensive MD simulations to investigate the structural properties of 
this model binary mixture as a function of a, varying this quantity in steps, 5a = 0.05, in 
the interval 0.2 < a < 0.5. For a = 0.25 we quantitatively related the interfacial properties 
with the corresponding phase diagram properties. In all the simulations we applied periodic 
boundary conditions along the x,y and z directions. At the lowest temperature, T* = 0.65, 
the simulations were initiated from a configuration where the molecules of type A and B form 
two contiguous FCC crystals. At higher temperatures we take as the initial configuration the 
final configuration of the previous temperature. The initial velocities of the molecules were 
assigned from a Maxwell-Boltzmann distribution. The equations of motion were integrated 
using a leap frog algorithm with a time step size, 6t* = 0.005. This corresponds to 1.1 x 10~^ 
nanoseconds in the scale of argon. At each time step iteration we monitor the temperature 
of the system, by means of the equipartition theorem, and rescale the linear momentum 
of the molecules to keep the temperature constant. This method of rescaling the linear 
momenta is known as the isokinetic thermostat. To check that this isokinetic thermostat 
produces meaningful results we have also carried out some MD simulations using the Nose- 
Hoover thermostat. In figure Q we show the kinetic energy distributions as the dynamics 
of the system evolves at, T* = 0.90, applying both methods. As one would expect both 
thermostats yield a Gaussian distribution of kinetic energies. The only difference is the width 
of the distributions. The isokinetic thermostat produced a distribution with a variance that 
is one order of magnitude smaller than the variance of the distribution obtained with the 
Nose- Hoover thermostat. As a further check of the isokinetic thermostat we also calculated 
and monitored some thermodynamic and surface quantities applying both thermostats. The 
results of this comparison test showed full consistency. Since MD simulations using the Nose- 
Hoover thermostat are computationally more demanding we used the isokinetic thermostat 
in all the simulations reported in this paper. Furthermore, to check the stability of the 
interfaces as well as the distribution of the species, we simulated the system for as long as 
55 ns in the scale of Argon. Thermodynamic quantities and interfacial properties of interest 
were measured by averaging over the last million of time-step iterations. To minimize 



correlations between measurements we calculated thermodynamic, structural and surface 
quantities every 50 time steps. We also investigated the role of finite size effects for the 
value of the miscibility parameter, a = 0.25. To this end we carried out MD simulations 
with three system sizes, N=1728, 4096 and 6144 molecules. We found that for all the 
quantities studied here, simulations with N = 4096 provided reliable results. Therefore, 
most of the simulations were carried out with = 4096 molecules. The discussion of 
these results will be presented, where appropriate, in the next section. On the other hand, 
interfacial properties are sensitive to the cross section area of the simulational box that is 
parallel to the interfaces, as discussed in previous MD simulations of the LV interface of a 



single LJ fluid |18|,| |19|. These authors arrived to the conclusion that a reliable value of the 



cross section area of the computational box should be at least (8cr)^. So, to be on the safe 
side, in the present simulations we have considered a computational box with a cross section 
area, Lx x Ly = (9(Taa)^- The length, Lz, of the simulational parallelepiped was adjusted 
such that the average density of the system laid somewhere inside the LV coexistence curve. 
In this way one readily gets the liquid-vapor phase coexistence. The average densities of 
the simulated systems were in the range 0.2 < p* < 0.4, where the reduced density defined 
as, p* = p(y\x- the following section we present, analyze and discuss the results of the 
thermodynamic and interfacial properties calculated from our MD simulations. 



IV. RESULTS AND DISCUSSION 



A. Structural properties and phase diagrams 

We performed extensive MD simulations for mixtures with 0.2 < a < 0.65 changing this 
parameter in steps of, 5a = 0.05. We studied the density profiles, p*{z), of the liquid- 
vapor coexistence in different temperature regions. From an analysis of p*{z), as function 
of T*, and, a, we estimated the Tw versus a phase diagram. In what follows we present 
some representative results for, p*{z), when a = 0.25 and 0.30. Then, for the particular 
value, a = 0.25, we present a quantitative analysis of the T* versus p*, and T* versus P* 
phase diagrams. We also give a brief description of the procedure we followed to locate the 
coexistence courve and the -A line- mixing- demixing line. 

Once the system reached equilibrium we calculated the structural properties of the system 
from simulations with N = 4096 molecules. In figure El we show the density profiles at the 



relatively low temperatures, T* = 0.65 and 0.75, when a = 0.25. In this temperature 
region the liquid-vapor equilibrium structure of the mixture consists of a liquid A-liquid B 
interface in coexistence with the vapor phase. As temperature increases and reaches the 
region 0.80 < T* < 1.25, however, this fluid phase structure rearranges in such a way that, 
the vapor phase spontaneously sets in between the Liquid A-Liquid B phases, as shown in 
figure 121 This structure remains stable during all the time of the simulation, about 55 ns 
in the scale of argon. This is a clear evidence that the vapor phase wets the liquid phases. 
In figure |3] we show similar results for a = 0.3. There one sees that at, T* = 0.82, the 
structure of the system is such that there is a liquid A-liquid B interface. Nonetheless, 
at higher temperatures, for instance, T* = 0.9, the vapor phase wets the liquid phases. 
This behavior of p*{z) for a = 0.3, suggest that T^(a = 0.25) < T^(a = 0.3). In fact, 
comparing the structures plotted in figures El and EJ one should note that the mixture with 
a = 0.25 already wets at, T* = 0.83. Following a systematic analysis of the structure of 
density profiles as function of temperatures for all the values of a we estimated the wetting 
transition temperatures T^{a). We found that increases monotonically as a function of 
a, whenever < a < 0.5. For higher values, 0.5 < a < 1.0, this wetting phenomenon does 
not occur. The results are summarized in the wetting phase diagram in figure El We believe 
that the reason for which the system wets below a = 0.5, and no longer does above this 
value is due to the equal size of the molecules of type A and B and because Eaa = ^bb- 

Now we will try to relate these density profiles structure and wetting phenomenon with 
the properties of the phase diagram of the mixture. To this end we will discuss in detail 
the properties of the corresponding phase diagram for a = 0.25. To begin it is important to 
remind that the present model binary mixture corresponds to the type III in the classification 
of Scott and Konynenburg ^20]. The phase diagram properties of this kind of mixtures has 



2l[, using a square well potential for the 



been quantitatively studied by Wilding et. al. 
intermolecular interactions. They showed that for a strong immiscible binary mixture a 
tricritical point exists. This means that the A line meets the LV coexistence curve at the 
critical point. In figure El we show the T* versus p* phase diagram of the mixture obtained 
from extensive MD simulations with a = 0.25. Since for a = and 1 our model reduces to 
single LJ fluids, we have also included in the same figure, for comparison, the phase diagram 
of a single LJ fluid. Both phase diagrams were calculated simulating a system with = 4096 
molecules and using a shifted intermolecular potential with a cutoff of 3cr. Therefore, the LV 
critical temperature of the single fluid became Tq"^ ~ 1.2. One should note that, the mixture 



critical density is higher than that of the single LJ fluid, and its critical temperature shifts 
upwards. This shift in p* and T* occurs due to the fact that the less miscible the mixture is, - 
smaller a-, the larger the temperature range of immiscibility. In addition, the A line appears 
to touch the LV critical point, and the critical point becomes tricritical. These result s ag ree 
well with those obtained recently for square well strong immiscible binary mixtures |2l|. 

To check how sensitive are the phase diagram properties of the mixture, for a = 0.25, 
to the number of molecules in the simulations, we also calculated the phase boundaries 
simulating a system with N = 1728 molecules. In fig. [7|we show the results of this finite size 
analysis. For a mixture with = 1728 the density of the liquid phase decreases while the 
density of the vapor phase increases. We also observe that at low temperatures the results 
are system size independent. However, as temperature increases, in particular, close to the 
critical point, there are differences in the phase boundaries obtained using = 1728 and 
4096 molecules. As one approaches the critical point it turns out more difficult to determine 
the coexistence densities with a system with A^ = 1728 molecules. This is so because the 
difference in the coexistence densities becomes smaller, and the number of particles in the 
system is not sufficient to give rise to bulk fluid phases. Nevertheless, for a system with 
A^ = 4096 this is not the case and we indeed obtained the liquid and vapor fluid phases. On 
the other hand, it is known that near to the critical point the fluctuations of the density 
are strong and the vapor and liquid densities are not well deflned. This fact complicates the 



location of the critical point. To circumvent this difficulty we proceed as in reference |22| and 
calculated the total density distributions in a system with A^ = 4096, at several temperatures 
around T^. The simulational box was divided in several slabs, parallel to the interface, of 
width between a and 3a . The density of particles was calculated in each slab every 50 time 
steps of the MD simulations. A block average histogram of the densities is obtained every 
5 X 10^ time steps. The resulting total density distribution, Pa + Pb, was calculated after 
averaging over 20 blocks. The result is presented in figure |H1. At temperatures slightly below 
Tc, we obtained density distributions that show two maxima. The low density maximum 
corresponds to the vapor phase and the higher maximum corresponds to the liquid density. 
Nonetheless, at temperatures above Tc, the density distributions show only one maximum. 
In FiglHlwe show the total density distributions at T* = 1.1, 1.15 and 1.2. The low density 
maximum is higher because the volume of the vapor phase is larger than the volume of the 
liquid phases. Thus, the vapor density appears with a higher frequency in the histograms. 
To locate of mixing- demixing transition temperatures, -A line- for a = 0.25 and A^ = 



4096, we followed a similar procedure as that described for the location of the LV critical 
point. However, in the determination of the A line, we only consider the density distribution 
of one of the species, pa or p^. The reason is that when the system is in the demixing 
region, the density distribution of one of the species shows two maxima. The low density 
peak corresponds to the poor fluid phase and the high density peak correspond to the rich 
fluid phase. Nonetheless, when the system is in the mixing region the fluid phases become 
homogeneous and therefore one should observe only one peak in the density distribution. 
Therefore, to locate the mixing-demixing points we looked at the transition from the two 
peak structure density distribution to one peak density distribution. This analysis was done 
as a function of the total density of the system and at three different temperatures. The 
results are shown in figure 01 

Another way of locating the A line is calculating the pressure versus temperature phase 
diagram, shown in figure EH for a system with N = 4096 molecules. The pressure was 
calculated as the average of the pressure tensor component perpendicular to the interface 
via the virial formula 23[. There we also included for comparison, the results for a single LJ 
fluid with the same number of molecules. We found that the LV phase boundary is located 
right at T* ^ 1.25. At higher temperatures there are two branches that were obtained by 
approaching the mixing-demixing boundary from both sides of the A line. These branches 
enclose a region that is narrower in size, as compared to the region obtained in the p* versus 
T* phase diagram. These results suggest that the calculation of the P* versus T* phase 
diagram gives more precise way to locate the mixing-demixing line. Again, we did find 
evidence that suggest that the A line meets the LV coexistence line at the LV critical point. 
These results give a strong support to the existence of a tricritical point j2l| . 



B. Interfacial properties 

Now we turn to the discussion of the interfacial properties and surface phenomenon of 
the mixture at the liquid-vapor phase coexistence. We carried out a quantitative analysis of 
these properties for a = 0.25. To evaluate the wetting transition temperature with precision 
we calculated the interfacial free energies as a function of temperature. To this purpose we 
use the well known formula, 

/•bulk2 , , 

7= / (Pn{z)-Pt{z))dz, (3) 

Jbulki ^ ^ 



where the integrations were carried out up to the middle of the corresponding bulk phases. 
The tangential and normal pressure profiles were calculated using the definition of the Irving- 
Kirkwood pressure tensor |23|. For a planar interface these pressure profiles are given by 



the formula 



24| 



P„(z) = piz)kBT (4) 
1 /^ 4<i(^»i) g|^ ^-^^ ^gpj 



2y4 \ ^ ^ Tij I Zij I V Zij J V 



Pt{z) = p{z)kBT (5) 



According to Young's rule the difference, 

A(T) = 27LV - 7LL, (6) 

must be zero at the WTT. So, the wetting by the vapor phase occurs when the free energy 
difference. A, becomes negative. This quantity is plotted as a function of reduced tempera- 
ture in the inset of figure ITT] for a = 0.25 and = 4096. A linear extrapolation of the data 
indicates that the wetting occurs at about T^(a = 0.25) = 0.80. Notice that due to the 
planar geometry of the interfaces there is no contact angle and the three surface tensions 
are independent and considered separately. To improve the accuracy of T^(a = 0.25), one 
needs to carry out even more demanding simulations. This is due to the fact that the in- 
terfacial tension always shows relatively large fluctuations. The situation complicates even 
more when the simulations are performed at temperatures very close to the WTT. A second 
independent way to estimate the WTT and to figure out the nature of the WT, we calculated 
the adsorption of molecules at the LL interfaces shown in figures El . This is done using the 
formula, 

/.bulkB 

r = / (p{z) - Phnlkjdz, (7) 

J bulkA ^ ^ 

The results of the calculations for a = 0.25 and A^ = 4096 and 6148 are plotted in figure IT^ 
One sees that r{T*,a = 0.25) decreases monotonically in the temperature range, 0.75 < 
T* < 0.78, and it consistently shows negative values. This is due to the inhomogeneity at 
the LL interface, since the density there, is much smaller than the density of the liquid bulk 
phases. Nonetheless, as T* —>■ 0.80 from below the adsorption jumps from a finite negative 
value up to minus infinite since the vapor wets the LL interface at T*{a = 0.25) ~ 0.80. Note 



that the WTT shifts shghtly towards higher temperatures as the number of molecules in the 
system increases from = 4096 up to = 6148. This jump or discontinuity is a strong 
indication that the WT is of first order. As expected the closer the temperature approach 
to T^(a) the stronger the fluctuations in r(T*,a). In figure IT^ the closest approach to the 
WTT, in reduced temperature, was ^ = 6.25 x 10^^. One would expect that the wetting 
transition is of first order even for other values of the miscibility parameter, < a < 0.5. 
This is so, since the behavior of p*(z, a) is similar to that of a = 0.25, whose surface 
properties were studied in detail. 



V. CONCLUSIONS 



We carried out extensive MD simulations to study the LV phase coexistence, the struc- 
tural properties and interfacial phenomena of a partially miscible symmetrical LJ binary 
mixture. By analyzing the density profiles as a function of temperature and miscibility pa- 
rameter we estimated the wetting phase diagram, versus a. The wetting of the vapor 
phase happens whenever < a < 0.5. We also found that T-^(a), monotonically increases 
up to a = 0.5. For other values of a, this wetting phenomenon does not occur. In addi- 
tion, we also studied quantitatively the T* versus p* and P* versus T* phase diagrams for 
a = 0.25. The results indicate that the former phase diagram shows a similar topology as 
that obtained for a square-well potential mixture estimated by means of mean-field theory 
and Monte Carlo simulations j21|. An analysis of the behavior of the adsorption of particles, 
at the LL interface, as a function of T*, for a = 0.25, led to the conclusion that the WT is of 
first order. These results should be valid for a family of mixtures of the type III, in the clas- 
sification of Scott and Konynemburg. The phase diagrams discussed here together with the 
wetting phenomenon are an explicit quantitative demonstration of the scenario suggested 
some time ago based on a microscopic expression for the Hamaker constant 2^ • To the best 
of our knowledge this is the first time that this wetting phenomenon is quantitatively studied 
by means of MD simulations in binary L J fluid mixtures in absence of a wall. We would 
like to point out that this wetting phenomena agrees well with some experimental studies 
carried out with fluid binary mixtures. Finally, the results reported in this paper provide 
with a more complete understanding of the surface phenomena in partially miscible fluid 
binary mixtures. 
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FIG. 1: Frequency versus reduced kinetic energy (Kinetic energy distribution) for a binary mixture 
with a = 0.25 and N = 4096 particles at T* = 0.9. The vertical axis should be multiplied by a 
factor of 10^. (a) Results using an isokinetic thermostat and, (b) results applying a Nose- Hoover 

thermostat. In each case the solid line represents the best fit to a Maxwell-Boltzmann distribution. 
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FIG. 2: Liquid-vapor-liquid reduced density profiles of the mixture for a. = 0.25 and N = 4096. 
Note that at these relatively low temperatures a liquid A-liquid B interface is formed. 




FIG. 3: Liquid-vapor-liquid reduced density profiles of the mixture for a. = 0.25 and N = 4096. 
At these higher temperatures the vapor phase spontaneously wets the liquid A-liquid-B phases. 




FIG. 4: Liquid-vapor-liquid reduced density profiles of the mixture for a = 0.30 
At these higher temperatures the vapor spontaneously wets the liquid A-liquid-B 



and N = 4096. 
phases. 
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FIG. 5: Reduced wetting transition temperature as function of a for an equimolar binary mixture. 
These results were obtained from MD simulations with a N = 4096 particles. The line is a guide 
to the eye. 
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FIG. 6: T* versus p* phase diagram for a partially miscible mixture with N = 4096 and a = 0.25. 
For comparison we have included the corresponding phase diagram of a single LJ fluid. 
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FIG. 7: T* versus p* phase diagram for a binary mixture with N = 1728 ( ♦ ) and = 4096 (0) 
particles and a = 0.25. 




FIG. 8: Frequency versus total reduced density (total density distribution) of the mixture with 
a = 0.25 and N = 4096 at T* = 1.1, 1.15, and 1.2. The vertical axis should be multiplied by a 
factor of 10^. These temperatures are slightly higher than the tricritical point temperature. 
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FIG. 9: Frequency versus reduced density (density distribution) for one of the fluid phases (fluid 
A) with a = 0.25 and N = 4096. The vertical axis should be multiplied by a factor of 10^. The 
left columns, show the density distributions when the A-line is approached form the mixing region 
at T* = 1.3, 1.4, and 1.5. Note that they show a one peak structure. However, if the A-Iine is 
approached from the demixing region at T* = 1.3, 1.4, and 1.5 the density distributions, right 
columns, show a two peak structure. 
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FIG. 10: Reduced pressure versus reduced temperature phase diagram for a = 0.25. 




FIG. 11: Reduced interfacial tension as function of reduced temperature for both, LL and LV 
interfaces. In the inset we plot the difference A = 27lv — 7ll as a function of temperature. These 
results are representative of a mixture with a = 0.25. 
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FIG. 12: Reduced adsorption as function of reduced temperature at the LL interface for two 
systems, N=4096 (A), N=6144 (y) and a = 0.25. The two points that are closer to the wetting 
transition temperature show relatively large error bars. This is due to the increase of fluctuations 
of the interface width close to T-^{a). In the present case, the vapor wets the LL interface at, 
T^{a) = 0.80, and the adsorption jumps from a finite negative value up to minus infinite. Dashed 
lines are a guide to the eye. 



